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ABSTRACT 

The first astrophysical objects shaped the cosmic environment by reionizing and heat¬ 
ing the intergalactic medium (IGM). In particular, X-rays are very efficient at heating 
the IGM before it became completely ionized, an effect that can be measured through 
the redshifted 21 cm line of neutral hydrogen. High-mass X-ray binaries (HMXBs), 
known to be prolihc X-ray sources in star-forming galaxies at lower redshifts, are prime 
candidates for driving the thermal evolution of the IGM at redshifts z ^ 20. Despite 
their importance, the formation efficiency of HMXBs from the first stellar populations 
is not well understood—as such, their collective X-ray emission and the subsequent 
imprint on the 21 cm signature are usually evaluated using free parameters. Using 
A-body simulations, we estimate the rate of HMXB formation via mutual gravita¬ 
tional interactions of nascent, small groups of the first stars (Population HI stars). We 
run two sets of calculations: one in which stars form in small groups of five in nearly 
Keplerian initial orbits, and another in which two such groups collide (an expected 
outcome of mergers of host protogalaxies). We find that HMXBs form at a rate of one 
per ^ 10^ Mq in newly born stars, and that they emit with a power of ^ 10^^ erg s“^ 
in the 2—10 keV band per solar mass per year of star formation. This value is a 
factor ^ 10^ larger than what is observed in star forming galaxies at lower redshifts; 
the X-ray production from early HMXBs would have been even more copious, if they 
also formed in situ or via migration in protostellar disks. Gombining our results with 
earlier studies suggests that early HMXBs were highly effective at heating the IGM 
and leaving a strong 21 cm signature. We discuss broader implications of our results, 
such as the rate of long gamma-ray bursts from Population HI stars and the direct 
collapse channel for massive black hole formation. 

Key words: cosmology: theory - early universe - cosmology: dark ages, reionization, 
first stars - stars: Population HI - intergalactic medium - X-rays: binaries - stars: 
kinematics and dynamics 


1 INTRODUCTION 


A major outstanding goal in cosmology is to piece together 
the history of the Universe between Cosmic Dawn, the emer¬ 
gence of the first stars and galaxies, and the end of reioniza¬ 
tion, when the radiation from these objects had ionized the 
intergalactic medium (IGM). Advances in numerical tech¬ 
niques, combined with exquisite measurements of the “ini¬ 
tial” conditions (at a redshift z « 1000; Hinshaw et al.|2013 
Planck Collaboration et al. 20151, have led to remarkable 


simulations (e.g. Abel, Bryan &: Norman|2002[ Turk, Abel Sz\ 
|Q’Shea|2009[ [Stacy, Greif fc Bromm|2010||Greif et al.|2011[ 
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[Bromm fc Yoshida|201l| ) of the conditions leading up to the 
former milestone, occurring at z > 30, when the Universe 
was « 100 Myr old. However, reconstructing the subsequent 
several hundred Myr of cosmic history has proved far more 
challenging, due to the difficulties in reliably modeling the 
numerous forms of feedback from the first astrophysical ob- 


jects (e.g. 

Springel, Di Matteo & Hernquist||2005t |Stinson| 

et al.|2006 

Sijacki et al. 120071. 


In particular. X-rays from the first galaxies can act as 
a powerful source of feedback (e.g. [Venkatesan, Giroux fcj 


Shull|[200l| Machacek, Bryan fc Abel||2003 l that exerts in¬ 

fluence over a wide range of distance scales. Because hard 
X-rays (energies > 1 keV) have mean free paths compara¬ 
ble to the Hubble horizon, they can isotropically heat and 
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partially reionize the early IGM (e.g. 

Oh|2001 

Venkatesan, 

Tumlinson & Shull|2003| [Ricotti & Ostriker|2004| [Pritchard 

& Furlanetto 

2007 

1. In fact, they are expected to be the dom- 


2010 Treister et al. 12012 1; because of the uncertainty in the 
fraction of X-ray photons that is released into the IGM as 
opposed to being trapped inside the accretion flow or re¬ 


inant agent in heating the IGM. Such heating may suppress processed into the infrared (e.g. Madau, Haardt & Dotti 


star formation | Ripamonti, Mapelli & Zaroubi 20081 and 


massive black hole (BH) growth (Tanaka, Perna & Haiman 


20121 inside low-mass dark matter haloes by raising the 
Jeans and filtering masses of the IGM ( Gnedin]2000 Naoz & 


Barkana||2007 1. On galactic and circum-galactic scales, soft 
X-rays (~ 0.1 — 1 keV) can affect the formation of stars and 
possibly massive black holes (BHs) by promoting the forma¬ 
tion of molecular hydrogen via electron-catalyzed reactions 
(e.g. Haiman, Rees & Loeb 1996| Kuhlen & Madau 2005 


Latif et al.||2015 Inayoshi fc Tanaka||2015 1. 

In addition to their suspected roles in early galaxy evo¬ 
lution, X-rays are important also because they can leave an 
observable signature that can be exploited to probe the cos¬ 
mological epoch in question ( [Pritchard &: Loeb|2008| ). Their 
thermal impact on the early Universe should be measurable 
through the redshifted 21 cm transition line of neutral hy¬ 
drogen, which is observed in emission or absorption depend¬ 
ing on the relative temperature of the IGM with respect to 
that of the cosmic microwave background (CMB). Several 
studies have investigated how forthcoming observations of 
the sky-average amplitude and power spectrum of the relic 
21 cm line (e.g. Bowman, Rogers fc Hewitt|2008| Burns et al. 


2012 van Haarlem et al.|2013 Voytek et al.||2014| ) could be 

used to constrain the astrophysical agent (or agents) respon¬ 
sible for heating the early IGM. 

There are sound reasons to expect that the first galax¬ 
ies produced X-rays in abundance, and rapidly heated the 
IGM. There are two dominant X-ray sources in present-day 
galaxies—both powered by gas accretion onto BHs, and both 
plausibly prominent shortly after Gosmic Dawn: gas feed¬ 
ing massive BHs shining as active galactic nuclei (AGN), 
and X-ray binaries, powered by a stellar-mass BH gradually 
cannibalizing a companion star. Estimates of the mass accu¬ 
mulated by nuclear BHs prior to z ~ 6 ([Shankar, Weinbergj 


& Miralda-Escude 2009 Salvaterra et al. 2012aI, the exis¬ 


tence of very massive BHs at a ^ 6 (e.g. 


Fan et al. 


well as the observed (e.g. Shen et al. 20071 and theoretr 
cally expected (Shankar, Weinberg fc Miralda-Escude|20TO 


2001 


McGreer et al.||2006| |Willott et 

al.[[ 20 UV| [ 2 UUy( 

u 

0 

u 

0 

et al.[[2011 

Venemans et al.[[2013 

Banados et al. 

2014 

), as 


Tanaka 20141 increase in their duty cycles toward higher red- 


shifts, all hint that X-ray AGN may have been much more 
common during this epoch. Likewise, high-mass X-ray bina¬ 
ries (HMXBs) dominate the X-ray emission of star-forming 
galaxies; the low metallicity and rapid baryonic mass accre¬ 
tion of the earliest galaxies both lend credence to the notion 
that they were rife with HMXBs. Theoretical models sug¬ 
gest that either type of X-ray source could heat the IGM to 
above the CMB temperature as early as « ~ 30, and that 
this transition should be measurable by the planned 21 cm 
experiments. 

At present, there are too many theoretical uncertain¬ 
ties to determine from the future data which type of X- 
ray source—AGN or HMXBs—was responsible for driving 
the thermal evolution of the 2 : < 30 IGM. Modeling the 
early AGN X-ray emission is particularly difficult, because 
the conditions for triggering AGN activity are not fully un- 


2014 Pacucci et al. 20151; and because the epoch, initial 
masses, and birthplaces of the massive BH “seeds” are not 


yet constrained by observations (e.g. Volonteri 2010 Haiman 
2013 Tanaka fc Li|2014 1 . Similarly, studies usually estimate 


the X-ray contribution from early HMXBs by simply infer¬ 
ring empirical relations between X-ray luminosity and star 
formation rate (SFR) in local galaxies (and modeling the 
SFR using a semi-analytic cosmological model), or combin¬ 
ing such relations with one or more free parameters (e.g. 
Mirabel et al.|2011| Tanaka, O’Leary fc Perna][2015 K 


The goal of this study is to alleviate the uncertainties in 
the formation rate and X-ray output of HMXBs in the early 
Universe, by using N-body simulations of nascent groups 
of the first (Population HI, henceforth Pop III) stars. We 
choose the properties of the star groups in our simulations to 
reflect those found in hydrodynamical simulations of Pop HI 
star formation at 2 ~ 20 (Greif et al.|2012 Stacy & Bromm 


20131. We follow the formation and dynamical evolution of 


compact binaries over thousands to millions of years, includ¬ 
ing the effects of the background gravitational potential and 
dynamical friction. This allows us to compute the fraction of 
Pop III stars that form stable, compact binaries, and even¬ 
tually undergo an X-ray bright phase. The end result is an 
estimate of the formation rate of HMXBs in the first pro¬ 
togalaxies, as well as the amount of X-rays they generate 
per unit star formation. To our knowledge, this is the first 
published estimate of this type. 

Our simulations predict a binary formation rate which 
is similar to what is observationally inferred in present- 
day galaxies. However, we derive a HMXB energy output 
(normalized to the star formation rate) that is a factor 
~ 10— 150 higher than in present-day star-forming galaxies, 
if the HMXB duty cycle is similar to the one in the local 
Universe. We find that the X-ray output does not change 
significantly within the wide variety of simulation setups 
considered—such as different orientations for collisions be¬ 
tween star groups, and ambient gas density—and submit 
that this is a robust estimate. 

The findings of this study can be used as model in¬ 
puts in estimating the 21 cm global signature and power 
spectrum, but have wider applications. As stated above, the 
X-ray output of the first galaxies are also of interest for 
studying feedback on smaller scales, such as subsequent star 
formation and massive BH formation. 

Our work is also relevant for predicting the rates of long- 
duration gamma-ray bursts (LGRBs) from Pop III stars. 
LGRBs are important probes that can shed light on the Uni¬ 


verse out to 2 > 10 (e.g. Toma, Sakamoto & Meszaros 20111. 
According to the collapsar model ([MacFadyen fc Woosley[ 


19991, progenitors of LGRBs require rapid rotation of the 


derstood even at low redshifts (e.g. Hopkins & Quataert 


He core and removal of the H envelope. Both criteria are sat¬ 
isfied by Pop HI HMXBs, and it is plausible that massive 
Pop HI stars in binary systems are dominant LGRB pro¬ 
genitors in the early Universe ( Bromm fc Loeb][2006 l. Our 
results on HMXB formation rates can therefore be used to 
predict and interpret observations of high-redshift LGRBs. 

The paper is organized as follows. We start in ©by 
discussing the problem to be solved—beginning with the 
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equations of motion, followed by the description of our N- 
body code, our choices for the initial conditions, and how 
the data is interpreted for HMXB formation. We present 
our results in In SQ we discuss the implication of our 
work for the X-ray output of the first galaxies, as well as 
for other topics such as LGRBs and SMBH formation. We 
conclude with a summary of our findings in ii 


2 STELLAR DYNAMICS 

Here, we provide an overview of our simulations—namely: 
the equations of motion that are solved to simulate the dy¬ 
namical evolution of the star groups; the numerical scheme 
we use to solve the equations; the different types of ini¬ 
tial conditions we adopted, as well as the reasoning behind 
our choices; and finally, how the results are interpreted for 
HMXB formation. 


2.1 The equations of motion 


Our A-body code computes the motion of N objects with 
(generally different) masses rui, moving under their mu¬ 
tual gravitational influence, a dissipative dynamical friction 
force, and a background gravitational potential. We numer¬ 
ically integrate the equations of motion 

<f ^ ^ ^ ^ ,, , 

^^2 ^bg,z- (1) 

The first term on the right-hand side of equation Q is 
the specific force due to Newtonian gravity. 


dsA — 


G ruj 


d S{rij) 

d ra 



( 2 ) 


where G is the gravitational constant, f) is the displacement 
of the ith star from the center of the host dark matter halo, 

and Tij = \fi — fj|. _ 

We adopt the Plummer softening kernel S{rij) (e.g. |Bin^ 
ney fc Tremaine|1987 1, 


S{rij) = - 


rfj + e^ 


(3) 


where we take e = 7?©. 

The second term on the right-hand side of equation Q , 
Odf.i, is the specific drag force due to dynamical friction. For 
collisionless systems, the standard Chandrasekhar formula 
for dynamical friction is ( Binney fc Tremaine|198T l, 

Si = -Itt InA f(Xi) —^ pifi) Vi, (4) 


where 


fiXi) = erf(W) - ^ W exp {-Xf) 


(5) 


Vi is the speed of the ith star with respect to the back¬ 
ground, Xi = Vi/{'/2av), a is the velocity dispersion, InA 
is the Coulomb logarithm and pifi) is the local gas density. 
We adopt the modified formula for gaseous medium 


used in Tanaka & Haiman (20091. This prescription incorpo¬ 


rates behaviors found in numerical simulations for subsonic 
and supersonic regimes ( Ostriker]|1999 Escala et al.||20d4|. 
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The specihc drag force vector always points opposite to the 
direction of motion, and is given by: 

^ ^ ( 6 ) 

with 



' 0.5 

InA 

erf 

^ Mi \ 

JlM, exp(-^)] 






0 s; Mi s; 0 . 8 ; 


1.5 

InA 

erf 


JlM^ exp(-:^)] 


0.8 ^ Af s; Afeq; 


+lnA 

M^ > Afeq. 

(7) 

Above, Mi = Vi/cs is the Mach number, and Ca is the sound 
speed. We use InA = 3.1 and the corresponding value of 
Meq « 1.5. 

In our simulations, the motion of the stars with respect 
to the background gas is supersonic. In this regime, the char¬ 
acteristic dynamical friction timescale, for a circular Keple- 
rian orbit of two bodies with for mi ';S> m 2 and V 2 Cs, 
is 


Eorh 1 / m\ 1 , . 

Pdf ^80YmiGp(i.,)i.3/2’ 

where Pdf is the frictional dissipation power (m 2 adt. 2 U 2 ), 
and Porb the orbital energy (Gmim 2 / 2 ri 2 ). 

The third and last term, flbg.i, is the specific force due to 
the background potential, which is dominated by gas. The 
background potential provides an additional inward force 
whose functional form depends on the density profile. For 
simplicity, here we use a constant density and explore dif¬ 
ferent values in our simulations. The force due to the back¬ 
ground potential is then 

flbg.i = -MGpvi, (9) 

where here Vi is the vector pointing from the center of the 
halo to the i-th star. 

The equation of motion is solved iteratively, with the 
positions, velocities and accelerations of each star updated 
at every time step. We describe our computational method 
below. 


2.2 Code Description 


We perform 3-dimensional, A-body simulations with 4th- 
order & 5-stage Runge-Kutta-Fehlberg methods (RKF45 
method, Erwin 19691 using adaptive time steps. The RKF45 


is a very precise and stable integration method among the 
large class of Runge-Kutta schemes, particularly by adapt¬ 
ing the Butcher tableau for Fehlberg’s 4(5) method. 

We solve equation 0- as described in the preceding 
text, updating the position and velocity components of the 
stars while treating the background density of the gas as a 
uniform and static distribution. 

To ensure numerical precision, our computational 
scheme varies the value of each subsequent time step an¬ 
alytically, so that numerical errors for each variable in the 
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simulation do not exceed 10“^® times the size of the vari¬ 
able. In some cases, however, this method leads to excessive 
computational effort for calculating relatively trivial inter¬ 
actions. For example, near the pericenter of hyperbolic or 
highly elliptical encounters, the time steps become increas¬ 
ingly small in a runaway fashion to compensate for the steep 
rise in acceleration and associated errors. In order to avoid 
such situations, we implement the following two numerical 
shortcuts to keep computation times tractable. 

The first shortcut is to use analytic approximations for 
very close 2-body encounters. This is justified in cases where 
pairs of stars are sufficiently close to each other and isolated 
from the other stars in the simulation, so that (i) the grav¬ 
itational pull from the other stars and the background po¬ 
tential are negligible compared to the mutual gravitational 
pull of the pair, and (ii) the orbital motion is supersonic 
and the dynamical friction force can be treated as a linear 
perturbation to the 2-body Keplerian problem. The code 
employs analytic approximations for any close stellar pairs 
that satisfy these conditions, and reverts to the RKF algo¬ 
rithm when the conditions stop being satisfied. 

We derived the following approximations for the semi¬ 
major axis and the eccentricity: 


a{t) 


o(fo) + P p 



-2/3 


(i - to) 


( 10 ) 




( 11 ) 


Above, p. is the reduced mass, a{to) and e(to) are respec¬ 
tively the semi-major axis and the eccentricity at t = to, P is 
a dimensionless constant, and w(q) is a function of the mass 
ratio that is symmetric about q = 1. Note that the variable 
a with subscript indicates the specific force, while without 
subscript it represents the semi-major axis. To derive the 
analytical expression for the time-evolution of the orbital 


distance in equation (101, we integrated the equation of the 
motion of a star under the influence of a dynamical fric¬ 


tion torque. Equation (111 then follows from the definition 


of orbital eccentricity in terms of orbital energy and angu¬ 
lar momentum. Note that as the radial distance decreases, 
the eccentricity increases. Here, p and w{q) are free param¬ 
eters and we use P = 0.035 and w{q) = q^ * -|- (l/qY'*. The 
stellar coordinates and velocities are recovered as functions 
of a and e using standard expressions for Keplerian orbits 
( Binney &: Tremaine|1987 |. 

Because we are dealing with systems with more than 
three bodies of different masses, stars frequently form hi¬ 
erarchical triple systems whose motions are affected by the 
Kozai mechanism ( |Kozai |T96^ . Since our analytic solutions 
above do not account for such changes in mutual inclina¬ 
tion, we limit the use of these solutions to situations where 
"Touter < TKozai, where Touter is the dynamical timescale for 
the outer pair in the hierarchical triple, and TKozai is the 
time scale for the Kozai mechanism. 


r)2 

^1 -^1,2 2 n1.5 

"?”Kozai ^ (i ^l,2j 

m2 Pi,3 


( 12 ) 


The subscript 1 above indicates the primary star of the inner 
compact binary along with the satellite star denoted by the 


subscript 3, while the hierarchical tertiary star is indexed by 
the subscript 2, and P is the orbital period. 

Our second shortcut for keeping the simulation run¬ 
times manageable is to set a minimum value for the time 
step. We choose a physically motivated value, 10~® x 
Tdyn, min 1 where Tdyn, min is the smallest value of the dynam¬ 
ical time between any two stars in the simulation (that are 
not being treated by the analytic shortcut above) at a given 
time step. This procedure is necessary when there are three 
or more stars interacting at small separations, in which case 
the analytic approximations above cannot be used. We note 
that such situations are rare compared to the places where 
the analytic shortcut is applicable. 


2.3 Determining HMXB formation 


It is assumed that a HMXB has formed if both of the fol¬ 
lowing criteria are satisfied: 

(i) One of two stars forming a binary turns into a black 
hole. In order to determine which stars turn into black holes, 
we need to compare the typical lifetime (rufe) of a massive 
star with the time (trun) in the simulation (taken to coincide 
with the time at which stars are born). If rufe > trun, the 
star is marked as a black hole in the simulation. For main- 
sequence stars, it is possible to estimate this lifetime from 
the mass-luminosity relation, that is Estav ~ m and Tatar ~ 
m^, leading to rufe ~ m}~^ with 2 < p < 3. However, due to 
the uncertain value of p for Pop HI stars, we rather prefer 
to use here the nuclear time scale of Pop HI stars estimated 


by Schaerer (20021 and Marigo, Chiosi & Kudritzki (20031. 


They calculate the H-burning nuclear time scale for these 
stars with a stellar evolution code. We assume that a BH 
forms if the stellar mass is greater than 8 M©. 

(ii) The two stars are close enough so that the accretion 
occurs through Roche-lobe overflow (RLOF). This criterion 
is simply written as Rstav S? Rrl, where Rstav is the stel¬ 
lar radius, and Rrl is the Roche-lobe radius of the most 
massive star calculated from the center of the star to the 
inner Lagrange point. An approximate analytic formula to 
the Roche-lobe radius of star 1 for a wide range of the mass 
ratio is ( |Eggleton| 1983P , 

^RL,! _ _ 0.49g^'^^ _ , , 

r 0.6g^/® -I- ln(l -|- g^/^) ’ 

where r is the orbital distance and g = — is the mass ratio. 
Although the above expression was derived in a study of 


binaries in circular orbits, Regos, Bailey & Mardling (20051 


found that the Roche lobe radius does not differ very much 
for eccentric binaries. We therefore apply equation ( |13| l to 
all binaries found in our simulation. 

If a binary satisfies these two criteria, the binary is 
marked as a HMXB. Since the two criteria are independently 
checked at every time step, which criterion is satisfied first 
is not important. 


2.4 Setup and Initial Conditions 

We design our simulations with the 20 < z < 30 Universe 
in mind. This is the redshift range where the dark matter 
haloes reach virial temperatures ~ 1000 — 2000 K, the ex¬ 
pected condition for Pop HI formation, at the highest rates. 
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The lifetimes of the stars are several Myr, which is com¬ 
parable to the timescales on which the host haloes undergo 
mergers with other Pop III forming haloes ( Tanaka|2014 |. 

The initial positions of the stars are generated quasi- 
randomly via a Monte Carlo realization as follows, moti¬ 
vated by the assumption that they formed inside a shared 
Keplerian gas disk at the center of the host dark matter halo. 
Their initial radial positions in the disk plane are chosen ran¬ 
domly from a uniform distribution whose amplitude depends 
on the characteristic size scale of the star-forming cloud (see 
below). Their azimuthal distribution is chosen to be nearly 
uniformly distributed, so that the azimuthal angular posi¬ 
tion of the n-th star is given by 360° (n/A'’) ± 5°, where N is 
the total number of stars in the group. This choice is made 
to minimize the radial gravitational pull between the stars, 
so that they do not immediately fly apart. Naturally, the az¬ 
imuthal positions do not remain evenly spaced, but become 
mixed quickly. Finally, we allow the stars to be displaced 
out of the plane of the disk. Their positions perpendicular 
to the plane are chosen randomly, so that their vertical dis¬ 
placement out of the disk plane is no more than 5% above 
or below their initial radial displacement from their shared 
center of mass. The initial velocities of the stars are assigned 
to be circular, parallel to the disk plane, and Keplerian at 
the instant the simulation begins. 

For this study, we investigated two different size scales 
of star forming regions. 


(i) Large scale: Following the results from Stacy & 
Bromm (20131, where Pop III protostars are formed inside 


star-forming regions with a size of a few thousand AU, we 
have performed simulations where the stars are placed inside 
a region of size ~ 2000 AU. 

We explore two different values of uniform, constant gas 
density, 10 ® cm~® (denoted ng hereafter) and lO'^ cm“® (n 4 ). 

For gas of primordial composition, the molecular weight 
H can vary between 0.6 and 1 . 2 , depending on the ionization 
fraction. We adopt ^ = 1 for simplicity; this choice does 
not qualitatively affect our results, since our values for n are 
selected arbitrarily with the goal of exploring the qualitative 
dependence on n. 

Each simulation is run for 5 Myr. The masses of the stars 
follow the initial mass function (IMF) they provide with a = 
0.17 (1^ = M-“), = 140 Mo and = 0.1 Mq. 

For these simulations, we consider N = 5 stars. 

(ii) Small scale: In the simulations by |Greif et al.| ( [2012[ ), 
multiple protostars formed several AU apart from each 
other. We therefore run a second set of simulations, where 
stars form within a 10 AU radius. Since [Greif et al.| ( |2012[ ) do 
not provide a slope for the IMF, we use a = 0.17 as above. 
That study found that more than half of the mass accreted 
during the protostar phase goes to the most massive proto¬ 
star in the group. 

To mimic this behavior, we first generate a star with 
Afmax = 200 Mq, then generate each subsequent star with 
Afniax = 200 Mq— [the sum of the masses of the previously 
generated stars]. 

Just as with the large-scale case, we run simulations with 
number densities ng and 714 . Because this case is the more 
relevant one for forming HMXBs, we explore several different 
configurations: cases with a star group of N = 5 stars, a star 
group of N = 10 stars, and collisions between two groups 


Large scale 

m 

714 

number of runs 

20 

23 

P{Bi 2 ) [P(Bi 2 -f B 13 )] 

0.60(0.90) 

0.56(0.74) 

(<lt=5Myr)[AU] 

270 

340 


Table 1. Summary of the large scale calculations for n = 
10 ® cm“^ (ug) and n = 10 ^ cm“® ( 714 ). P(Bi 2 ) denotes the frac¬ 
tion of the simulations in which the most compact binary consists 
of the two most massive stars. We also list P(Bi 2 -I-B 13 ), the frac¬ 
tion where the most compact binary consists of the most massive 
star paired with either the second or third most massive star. Also 
shown is the average semi-major axis of the most compact binary 
at t = 5 Myr. 


of N = 5 stars. The last case is motivated by the fact that 
Pop Ill-forming minihaloes undergo frequent mergers, which 
suggests that the nascent star groups themselves undergo 
close encounters. 


3 RESULTS 

Here, we summarize the findings of our Af-body simulations, 
focusing in particular on the properties of the most compact 
binaries found for each set of runs. 


3.1 Large scale 


We ran 20 simulations for 77 g and 23 for 714 . 

For ng, in 12 out of 20 runs, the two most massive stars 
(Si and S 2 , where Si is the i-th most massive star) form the 
most compact binary—we denote such a binary with the 
notation B 12 . Similarly, binaries of type B 13 (i.e. made up 
of the most massive star Si and the third most massive star 
S 3 ) form the most compact binary in 6 of the runs. 

The left panel of Figure shows a sample set of stellar 
trajectories for one of the large scale Tig runs. Throughout 
the simulation, the compact binary tends to remain near the 
center of the halo, as less massive stars repeatedly undergo 
3-body interactions with the binary. The fact that they do 
not stray far from the center of the halo is a combined effect 
of the background potential and the gravitational potential 
of the massive binary. It is these 3-body encounters that 
cause the most compact binary in the simulations to end up 
as type B 12 or B 13 . In two of the runs, the most compact 
binaries after 5 Myr consist of less massive stars. 

The average value of the semimajor axis after t = 5 Myr 
in the large scale. Tig runs is Myr) = 270 AU, and the 
minimum value of a across the 20 runs is 60 AU. These val¬ 
ues for a are much larger than that necessary for RLOF to 
take place, ~ 0.07 AU. Their characteristic dynamical fric¬ 
tion time scales are roughly 10 ^® yrs, which is much longer 


than their lifetimes (see equation 13 when Prl = Rstar and 
equation|^. We therefore conclude that in systems of A^ = 5 
stars and n ~ 10 ® cm“®, with the stars initially separated 
at hundreds of AU, HMXBs are unlikely to form. 

For the 714 case, the fraction of runs where the most 
compact binary after 5 Myr is type B 12 or B 13 are similar 
to the Tig simulations: 13/23 for B 12 , and 17/23 for B 12 or 
B 13 . 

However, the dynamical evolution is quite different, as 
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Figure 1. Left-. Sample trajectories for large scale calculations with the high number density n = 10® cm“® (case ng). As shown in the 
figure, stars tend to stay near the center of the halo and their overall motions are oblate-spheroidal in shape. Right. Sample trajectories 
for the low number density n = 10^ cm“® (rn). Even though stars including binary systems remain within a certain distance range, 
they are not as close as the stars in the higher density calculation, leading to less frequent three-body interactions. Furthermore, one can 
notice that a few stars (blue line and purple line) are kicked off. 



t [yrl 

Figure 2. Evolution of the semi-major axis of a typical binary 
in a system. The sharp variations are due to stellar scatterings, 
which mostly result in a hardening of the binary. After the multi¬ 
ple system gets stabilized and isolated (which happens after Ri 10 
years), the decrease rate of the semi-major axis depends on dy¬ 
namical friction alone. 

can be seen by the right panel of Figure The lower gas 
densities lead to weaker forces due to dynamical friction 
and background potential, and stars (especially less mas¬ 
sive ones) tend to be scattered farther from their initial 
position, as far as > 1 pc. This also results in less fre¬ 
quent three-body interactions in general, and may explain 
the larger average value of the semi-major axis after 5 Myr, 
(at= 5 Myr) = 340 AU. The formation of HMXBs appears even 
more intractable for the n 4 case. 

Sample stellar trajectories for the ne and n 4 large scale 
cases are shown in Figure and a summary of the results 
is given in Table 

3.2 Small scale 

For the small scale case, we performed 86 runs of 5-body 
simulations for each of the number density values ne and 


Small scale 

UQ 

724 

runs 

86 

86 

(a<=500yr} [AU] 

1.37 

1.42 

T-dfiyr] 

~ 10 ^® 

~ IQi® 

Companion stars 

3rd massive star. 

, 11 ~ 12 M© 

FhmXBc 

0.070 

0.070 

Fhmxb [10“"^ M“^] 

4.6 

4.6 


Table 2. Summary of results for simulations of 5-body groups 
forming on small scales, {at —5000 yr) indicates the semi-major 
axis at t = 5000 yr, while rjf represents the dynamical friction 
timescale required for at=5000 yr to shrink to apiL (see equation 
[^- FjjmxBc is the fraction of runs in which a HMXBc forms, and 
Fhmxb is the number of HMXBc formed across all simulations, 
normalized by the total mass of the stars in the simulations. 


n 4 . Due to the smaller initial separations of the stars, we 
run the simulations for a shorter amount of time. 

We run each simulation for a minimum of 1000 yr, but 
stop the run if a stable binary forms, and if no further signif¬ 
icant dynamical changes are observed. If such a binary does 
not form, we run the simulations to a maximum duration of 
5000 yr. 

In order to properly compare the results from the ne 
calculations with those from the n 4 calculations, we use the 
same initial conditions for each set of runs. 

In the ne calculations, the most common scenario is 
that Si always forms the most compact binary almost im¬ 
mediately, while stellar scatterings are most common during 
the first few years to about 40 years. Thereafter, a multiple 
system usually survives and stabilizes, while less massive 
stars are ejected. A difference between the large-scale and 
the small-scale scenarios is that, whereas most cases in the 
large-scale calculations end up with one binary and unbound 
single stars (in 15 out of 20 cases), in the small-scale runs a 
multiple system such as triple or quartet (rather than a sim¬ 
ple binary) forms (in 59 of 86 cases) Less massive stars take 
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Figure 3. Left: Pericenter distance - eccentricity distribution plot for riQ. Two HMXB candidates (HMXBc) have been produced. Right: 
The same distribution plot for 714 . The circled points indicate HMXBc whose mass transfer via RLOF may occur periodically due to 
their eccentric orbits. The elliptically-circled point indicates an HMXBc whose semimajor axis is smaller than uitl, so the mass transfer 
will be steady. There are a couple of binaries with pericenter distance smaller than Rrli but they are excluded because M 2 < 8 Mq. 


some energy from the multiple system and convert it into 
their kinetic energies while causing the binary to harden. 

As an example, Figure shows the evolution in the 
semi-major axis of a compact binary in one of the runs. One 
can easily notice the quick decrease in the semi-major axis 
during violent stellar interactions (before 10 years), which 
clearly indicates the hardening. After a couple of stars are 
cast away and the multiple system is stabilized and isolated 
(after about 10 years), dynamical friction plays the main 
role in the decrease of the semi-major axis; after this point, 
we do not expect further dramatic hardening of the binary. 
In our 86 runs, the last surviving binary is typically of type 
Bi 3 (pairing of the most massive and third most massive 
stars), with total mass of 130 M© and {Msfi = 11 M©. 
(at= 5 ooo yr) = 1-87 AU and the minimum semi-major axis 
is 0.2 AU. The corresponding dynamical friction time scale 
(equation]^ is ~ lO'^^ yr. The ejected stars have speeds of 
(10—100) Cs , and Cs ~ 4 km/s. Six of the HMXB candidates 
(HMXBc hereafter) have been formed in all runs (Phmxbc 
= 0.070). For later use, let us define Phmxb as the number 
of HMXBc formed across all simulations, normalized by the 
total mass of the stars. Then Phmxbc = 4.6 x lO”'^ ^ 0 ^ 
(This term will be used to estimate the X-ray luminosity 
and is one of the primary results of our study). 

The outcomes of the 724 simulations are quite simi¬ 
lar: a triple or higher multiple forms in 65 of 86 runs, and 
(at= 5 ooo yr) = 1.42 AU. Furthermore, due to the same num¬ 
ber of HMXBc, Phmxbc and Phmxbc are the same. One 
notable difference is that for n 4 , the dynamical friction time 
scale is longer by 2 orders of magnitude compared to ng, be¬ 
cause this quantity is inversely proportional to the number 
density. 

Interestingly, there are four runs (for each density value) 
in which the most compact binary is eccentric, and inside 
the requisite separation for RLOF at pericenter but outside 
it at apocenter. We consider only two of them as HMXBc 
and rule out the other two binaries since the mass of the 
more massive star is smaller than 8 Mq ( Heger et al.|2003 |. 
We present the distribution of eccentricity for pericenter dis¬ 
tance in Figurej^ In particular, the left panel shows the dis¬ 


tribution for Tie calculations and the right panel for 114. The 
circled point indicates HMXBc with distance at pericenter 
shorter than the corresponding Roche-Lobe radius. 

How accretion proceeds in a highly eccentric binary sys¬ 
tem under these conditions remains an unsettled issue to 
date, as studies have claimed that the orbital semimajor axis 
and eccentricity can either increase or decrease depending 


on the binary properties at pericenter (Sepinsky, Willems & 
|Kalogera|2007| and |Sepinsky et al.|2009[ ). If the RLOF does 
induce circularization, then accretion proceeds normally (i.e. 
steadily). However, if the RLOF instead increases the eccen¬ 
tricity of the system, then, whether accretion can proceed 
steadily rather than intermittently will depend on the rela¬ 
tive timescale between the disk lifetime Tdisk (on the order 
of the viscous timescale), and the orbital period of the bi¬ 
nary torb- We computed these timescales for all the eccen¬ 
tric binaries in our simulations (for which the Roche-Lobe 
radius straddles the pericenter and apocenter), and found 
that Tdisk > torb iu all cases but one. 

This implies that the fraction of binaries whose eccen¬ 
tricities cause intermittent RLOF is small, and that as a 
global average, RLOF is steady to a good approximation. 

The magnitude of the gas density considerably in¬ 
fluences the characteristics of the dynamical interactions. 
While on the one hand 40% of the simulations end up form¬ 
ing the same triples for both density values, on the other 
hand their trajectories and the center of mass movements 
differ significantly. Figure|^shows two sample trajectories of 
5 stars with low number density (left panel) and high num¬ 
ber density (right panel) after 1000 yr. They were given iden¬ 
tical initial conditions for the run, but their trajectories have 
developed differently. In Figureat 800 yr, even though the 
same stars form a triple and the same star is ejected for both 
number densities, their trajectories are clearly different. 

As can be seen in Figure which shows the change of 
forces per unit mass exerted on a typical kicked-off star, Udf 
and flbg are negligible compared to Ugr before about 5 years. 

Note that Odf and Ubg are synchronized at early times 
because the motions of the stars are close to the Keplerian 
motion, r ~ 1/u^. The star in Figure|^is ejected at t « 10 yr 
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Figure 4. Trajectories up to t = 2000 yr for 5-body simulations with identical initial conditions but different densities: n 4 (left panel) 
and ng (right panel). Despite having the same starting point, the trajectories of the 5 stars evolve quite distinctly in the two backgrounds 
due to the different magnitudes of dynamical friction. The two simulations for different bound systems: a binary for 714 (yellow line) and 
a quartet for ng (red-|-blue-|-yellow lines). 



Figure 5. Trajectories of two 5-body simulations with identical initial conditions but different ambient gas densities. The simulations 
have formed triples with the identical stars (yellow+purple lines; contrast with[^ at t = 1000 yr, despite differences in the trajectories 
of each star and the center of mass. 


for the high number density (and at 60 < t < 80 yr for the 
low nnmber density). This can be nnderstood from the fact 
that ttgr monotonically decreases and ogg increases (agg ~ r 
where r is the distance from the center of mass). At the 
same time, agf barely changes since the star moves in the 
supersonic regime and the speed decreases slowly. 

There is no noticeable difference in the overall results 
between the and ng runs—quantities such as (at^sooo yr), 
the dynamical friction time scale rgr, the total mass of the 
most compact binary (or which star forms the compact bi¬ 
nary with Si) as summarized in TableFor both number 
densities, the companion star of the binary is typically the 
third most massive star S3. 

To sum up, we find that Phmxbc ~ 7% of our simu¬ 
lations form HMXBc, regardless of the gas density value. 
Normalized to the total stellar mass in the simulations, the 
number of HMXBc formed per stellar mass is Phmxb ~ 
4.6 X 10""^ Mq\ 


3.3 10-body simulations 


We now explore several different configurations for the star 
group, and run several sets of simulations with 10 stars (in¬ 
stead of 5). These are: (1) 10-body version of the small scale 
calculation presented above; (2) head-on crash of two star 
groups containing 5 stars each; (3) a close encounter and 
subsequent inspiral and merger of two star groups contain¬ 
ing 5 stars each. The latter two scenarios are motivated 
by the fact that the merger timescales and mass accretion 
timescales of Pop III host haloes, as well as the lifetimes of 
the massive Pop III stars themselves, are of the same or¬ 
der, ~ 10 Myr. This suggests that merging haloes will be 
continuously forming new stars (perhaps Pop II instead of 
Pop III) as they merge with other haloes, and that close 
interactions and mergers of nascent star groups may be rel¬ 
atively common. We generate stellar masses in the same way 
as for the 5-body case, but with a larger value for the pa¬ 
rameter Mmax = 300 M©. We have run each simulation for 
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t [yr] 

Figure 6. Changes with time in the three forces per unit mass 
(gravity from other stars, dynamical friction and background 
gravity) acting on star 4 (the light blue color line in Figure]^. 
We run two simulations with identical initial conditions but dif¬ 
ferent ambient gas densities 714 and ng. During the early phases 
of the interaction, the gravitational force (ugr) dominates by sev¬ 
eral orders of magnitude over dynamical friction (auf) the 
background force (obg)- Note that agr,ne (black line) and agr,n 4 
(green line) synchronize at early times because the motion of the 
star is very close to Keplerian and r ~ l/u^. After the star is 
ejected at 10 yr for ng (later for 714 ), Ugr monotonically decreases 
and agg increases (agg ~ r where r is the distance from the center 
of the halo). During the same time, agf barely changes since the 
star moves in the supersonic regime and the speed decreases very 
slowly. 


Scenario 

1 

2 

3 

4A 

4B 

runs 

54 

30 

30 

30 

30 

(at=500 yr) [ AU] 

1.0 

1.1 

0.90 

1.8 

1.6 

(UBI^) [AU] 

1.6 

1.5 

1.8 

2.4 

2.0 

(a.e.t) [ AU] 

0.72 

0.15 

0.42 

1.2 

0.23 

FhmXBc 

0.33 

0.33 

0.27 

0.13 

0.27 

Fhmxbc [lO”'^ M“^] 

15 

11 

9.0 

4.2 

8.4 


Table 3. Summary of the results from 10-body simulations. We 
have considered five different scenarios. Scenario 1 simulated the 
10-body version of the 5-body calculation, i.e. an isolated star 
group. The other scenarios all involve collisions of two 5-body 
star groups. Scenario 2 is a head-on collision of two coplanar, co¬ 
rotating star groups. Scenario 3 also collided two groups of stars 
with co-rotating orbital planes, but with an impact parameter 
comparable to the sizes of the groups, which results in an inspiral 
and eventual merger. Scenarios 4A and 4B are similar to scenario 
3, except that the orbital planes of the colliding star groups had 
mutual inclinations of 45° and 135°, respectively. 

500 yr. The results of these simulations are summarized in 
Table We briefly discuss each one, as follows. 

3.3.1 Scenario 1: 10-body group in isolation 

We set up the simulations as in the 5-body calculations, but 
with 10 stars. 

We find that (at^sooyr) = 1.0 AU. Interestingly, there 
are 18 out of 54 cases in which a < orl. There is a large 
difference in scale between auia Orest, where auia is the 


semi-major axis of B12 (binary made up of the two most 
massive stars) and Orest is the semi-major axis of the binary 
stars other than Si and S2. This is a common feature of the 
10-body simulations: they often end up with triples whose 
inner binary is Brest while the outer binary is B12. Our sim¬ 
ulations yield (orest) = 0.72 AU and (0312) = 1-6 AU. Also 
note that, in 14 out of 18 HMXBc, the binary is Brest (i.e. 
it is not made up of the two most massive stars). 

Since the compact binaries in these simulations form 
quickly and we only follow them for 500 years, it is techni¬ 
cally possible that they will be disrupted before one of the 
stars turns into a BH. However, our simulations for the 5- 
body scenario showed that compact, quasi-steady binaries 
are unlikely to be disrupted, and for practical purposes we 
extrapolate this qualitative result to the 10-body case. 

We find that a HMXBc forms in a larger fraction of 
these simulations than in the 5-body case, Phmxbc = 0.33, 
for the obvious reason that there are more stars. Per unit 
stellar mass in the simulations, the number of HMXB can¬ 
didates is Phmxb = 1.5 X 10“^ ^0^’ which is a factor « 3 
higher than we found for the 5-body case. 

3.3.2 Scenario 2: Collision between two 5-star groups - 
head-on collision 

Two groups of 5 stars are set up with random initial condi¬ 
tions, in the same manner as for the previous simulations of 
5-body groups. The two groups are then arranged to collide 
head-on, as follows: they are placed at a separation of two 
to three times their sizes and their disks are aligned so that 
the mutual inclination is zero. The initial relative speed of 
the groups is roughly the speed of sound and the center of 
mass of one group is set to move directly toward the center 
of mass of the other group. 

We find that prior to colliding, each group forms a com¬ 
pact binary of type B12 (the most and second-most massive 
star). When the two groups collide, those two binaries that 
existed before the collision were broken and the two most 
massive stars of each group form a new compact binary with 
high chances. 

The average (04=500 yr) is 1.1 AU, but (obij) = 1-5 AU 
and (orest) = 0.15 AU, meaning that the most compact bina¬ 
ries are not formed from the most massive stars. The shorter 
average separations may be a result of a larger number of 
early 3-body scatterings, which act to harden the group as 
a whole. 

HMXBc form in 10 out of 30 runs, and they are not of 
type Bi 2. However, we find a rate of HMXB formation per 
stellar mass Phmxb = 1-1 x 10“® Alg^; this is higher than 
in the 5-body case and comparable to Scenario 1 above. 
Sample trajectories for one of the simulations of scenario 2 
are depicted in Figure 

3.3.3 Scenario 3: Collision between two 5-star groups - 
spirally merging case 

We have used the same input parameters for the two groups 
as in Scenario 2, except that we now set the impact param¬ 
eter to be of order the size of the group, whereas it was set 
to zero in Scenario 2. 

The groups are given opposite velocities of ~ Ca along 
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Figure 7. 10-body head-collision (scenario 2). In these sample trajectories, five stars in each group are in nearly Keplerian motion. The 
groups move towards each other with a relative velocity of ~ Cs (left panel). At t ^ 8 yr (middle panel), the two haloes begin to merge. 
During the merger, all ten stars undergo stellar scatterings (right panel). One triple (black-l-dark blue-|-red lines) has been formed and 
it is moving in -y direction. In these head-on collisions, it is likely for multiple systems that existed before merging to be broken, while 
a few new multiple systems are formed. 



Figure 8. 10-body spirally merging case (scenario 3). As with scenario 2, two groups of 5 stars are set on a collision cours. However, 
in this scenario the impact parameter is ~ 2 — 3 times the size of groups, whereas in scenario 2 it is set to be zero. The groups are 
approaching each other at the relative velocity of Ca (left panel). At t = 8 yr, the groups are about to merge (middle panel). After some 
time (right panel), two binaries are ejected (red and purple lines in the x-direction and black and blue lines in the j/-direction). In this 
sample case, the most compact binary is the one ejected in the y direction (black+blue lines) with a = 0.8 AU at f = 100 yr. A difference 
between this scenario and the head-on collision is that compact systems are more likely to survive the merger. 


the a;-direction, and are offset by a displacement along the 
y-direction that is ~ 2 — 3 times the typical size of the star 
group (~ 20 AU) so that they merge with a spiral motion. 

We find that (ot^soo yr) = 0.90 AU, (0812) = 1-8 AU 
and (arest) = 0.42 AU. 

The average separation lies between what we find in Sce¬ 
narios 1 and 2. This can be interpreted as being due to the 
fact that these simulations (in which the two groups merge 
gradually via inspiral) have more close 3-body interactions 
than in Scenario 1 (in which 10 stars in quasi-Keplerian or¬ 
bits evolve in isolation) but fewer such interactions than in 
Scenario 2 (in which the two groups merge head-on). 

HMXBs form in 8 out of 30 runs and, as with Scenarios 
1 and 2, none of the HMXBs are made up of the two most 
massive stars. We find a similar HMXB formation rate per 
stellar mass, Uhmxb = 9.0 x 10“^ ^0^- Sample trajecto¬ 
ries from one of the simulations for Scenario 3 are shown in 
Figure 

3 . 3.4 Scenarios 4^ and 4B: Collision between two 5-body 
groups - spirally merging case with inclinations of 
45 degrees and 135 degrees 

In these two scenarios, we again set two groups of five stars 
each on a collision course. The difference is that the orbital 
plane of one star group is tilted, so that the two stellar disks 


have a mutual inclination i. We set the inclination at i = 45° 
(nearly co-rotating) for scenario 4A, and i = 135° (nearly 
counter-rotating) for scenario 4B. The groups are initially 
placed at a separation of two to three times their sizes, and 
set in motion at the same speeds as for the inspiral case 
(scenario 3). 

We find that (at=5oo yr) = 1-8 AU, (ubij) = 2.4 AU, 
and (urest) = 1.2 AU for scenario 4A. In scenario 4B, we 
find (at=500 yr) = 1.6 AU, (obis) = 2.0 AU, and (urest) = 
0.23 AU. In 70% of the runs for both scenarios, the two most 
massive stars form the most compact binary. Stellar bina¬ 
ries end up with somewhat closer separations in the nearly 
counter-rotating case, due to the fact that the net angular 
momentum of the merged star group is smaller. Indeed, we 
find a total of four HMXBc across all the i = 45° simula¬ 
tions, and eight in a same number of i = 135° simulations. 

For the same reason, we find a larger fraction of HMXB 
candidates per stellar mass simulated in the nearly counter¬ 
rotating case (Fhmxb = 8.4 x 10“"* Mq) compared to the 
nearly co-rotating case (Fhmxb = 4.2 x 10“^ Mq). The 
overall formation rate of HMXB candidates is lower for both 
cases than the cases in which all the stellar orbits were nearly 
coplanar (Scenarios 1, 2 and 3), plausibly due to the addi¬ 
tional degree of freedom in the stellar orbits. Still, the value 
of Fhmxb is within a factor of a few for all of our simulations. 
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Figure 9. Trajectories for scenario 4A, a collision between two 5-body groups, each on quasi-Keplerian orbits, but with the orbital plane 
of the two groups tilted with respect to each other at an inclination i = 45°. The setup is the same as scenario 3, except for the mutual 
inclination of the orbital planes of the colliding star groups. The left panel shows the two groups on a collision course. At t = 4 yr, the 
halos are about to merge (middle panel). After t = 10 yr (right panel), a triple having the most compact binary are ejected [black and 
blue lines (inner binary) and brown line]. 



Figure 10. Trajectories for scenario 4B—the same as scenario 4A (Figure]^, but with the two groups nearly counter-rotating with 
respect to each other (f = 135°). At t = 4 yr, the groups are about to merge (middle panel) and after t = 10 yr (right panel), a quartet 
containing the most compact binary is ejected toward the upper right of the panel [black-|-blue lines (inner binary), brown and thick 
purple lines]. 


Sample trajectories from runs for scenarios 4A and 4B 
are depicted in Figure and Figure respectively. 


4 DISCUSSION 

4.1 Binary evolution and formation of HMXB 
candidates 

Our large-scale simulations show that the if Pop III stars 
form hundreds of AU apart, as in the simulations of [Stacy] 
& Bromm (20131, then the timescales required to make 


HMXBs via stellar scatterings are simply too long. On the 
other hand, if protostellar clouds fragment and form stars in 
close groups on scales of ^ 10 AU, as in Greif et al. (20111, 


then a small fraction of groups can form HMABs. We brietly 
discuss the dynamics of HMXB formation in our simulations, 
then move on to discuss the astrophysical implications of our 
Hndings. 

The simulations indicate that, as expected, scatterings 
play a major role in making a compact binary. The back¬ 
ground potential and the dynamical friction play a sec¬ 
ondary role, by allowing the most compact binary (or triple) 
to remain near the center of mass of the halo and for other 
stars to return and scatter again and again. 

We Hnd that on average, the number of HMXB candi¬ 
dates formed per stellar mass, Fhmxb, is a function of the 
number and orientation of close 3-body encounters. Our re¬ 


sults indicate that Fhmxb may be somewhat higher in con- 
hgurations that result in fewer ejections of stars, and if the 
interactions are coplanar. While we are able to interpret this, 
as well as trends in the average separation between stellar 
pairs, in terms of the initial kinematic setup of the various 


scenarios simulated (see [3.3 above), the value of Fhmxb 
does not vary by more than a factor ~ 3. We interpret this 
lack of a significant variation in Fhmxb, for such a diverse 
set of initial conditions and ambient gas densities, to mean 
that our values are not far from the one that results from 
similar stellar encounters in nature. 


4.2 The effect of migration on the formation of 
HMXBs 

Another way in which a nascent stellar group could harden is 
migration through a gaseous disk. The migration could occur 


as the protostars form—Greif et al. (20121 found significant 


accretion from the protostellar disk onto the most massive 
protostar, and did not follow the evolution of the system 
beyond this stage. (It could also occur to a lesser degree in 
a vestigial gas disk, after the stars are in place.) 

We evaluate the possible role of disk migration on the 
separation of Pop HI stars by considering a steady, geometri¬ 


cally thin disk with an a viscosity (Shakura & Sunyaev 1973 
see also [Frank, King fc Raiiie||2002[ ). We adopt a disk with 
a = 0.01 and an accretion rate m ~ 10“^ Mq yr“^, follow- 
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ing Tan & McKee (20041 and Tan & Blackman (20041, who 


considered the structure of accretion disks around Pop III 
stars at high redshifts. 

We estimate the migration timescale Tmig following|Syer| 


& Clarke (19951. We take a binary system with primary 
mass Ml = 120 Mq and secondary mass M 2 = 11 M©, 
based on the mean values found across our simulations. For 
these masses and disk parameters, the secondary is able to 
clear a gap around its orbital path (e.g. Syer fc Clarke|19^ 


|Seager|[2010[ [Lubow fc Ida|[^010[ ) , by satisfying both of the 
following two conditions: 




1/5 


(14) 


(15) 



10 
a [AU] 


15 


20 


where R is the distance of the secondary from the pri¬ 
mary, H is the scale height of the disk at that location, 
and q — M 2 /M 1 ~ 0.1 is the binary’s mass ratio. If con¬ 
dition (1) is violated, the gap will be closed by the radial 
pressure gradient. Condition (2) relates the gap width and 
the Roche radius of the secondary; this ensures that the sec¬ 
ondary acts to transfer orbital angular momentum through 
the disk, rather than accreting mass via RLOF. We find that 
inside 7? ~ 15 AU, for the disk properties stated above, H/R 
is smaller than the right-hand side of equation (141 by a fac¬ 
tor ^ 80, and smaller than the right-hand side of equation 
(151 by a factor > 8, indicating that the secondary is easily 
able to open a gap. 

The migration timescale of the secondary depends on 
the dimensionless parameter (|Syer &: Clarke|1995|) 


B = 


47rEoR^ 

M2 


(16) 


where Eo is the local surface density of a steady-state disk 
around the primary, in the absence of perturbations by the 
secondary. For B > 1, the gas in the disk is able to dy¬ 
namically dominate over the gravitational influence of the 
secondary, and the secondary is pushed inward on the vis¬ 
cous diffusion timescale of the disk. 


"^mig,0 



(17) 


For B < 1, Syer & Clarke (19951 found that the migration 
timescale is longer. 


Tmig,! ~ ^7/17 ■ (18) 

In Figure [ni„ we plot the local disk mass near the 
secondary, 47rEoi?^, alongside the typical secondary mass, 
11 Mq (the numerator and denominator, respectively, for 
the ratio B). For i? < 15 AU, B < 1 and the migration 
timescale is expected to slow. We plot the two migration 
timescales Tmig.o and Tmig,i in Figure [T^ 

Both of these timescales are shorter than both the typ¬ 
ical lifetimes of protostellar disks, as well as of the stars 
themselves. This suggests that disk migration could, in prin¬ 
ciple, lead to initial stellar separations smaller than what 
we have assumed, making the formation of HMXBs via stel¬ 
lar scatterings more favorable. On the other hand, radiative 
feedback from the stars could blow away the disk before sig- 
nihcant migration can occur. We submit that the HMXB 


Figure 11. A comparison between the mass of a typical sec¬ 
ondary of a binary in our simulations (M 2 = 11 Mq, dashed 
line) and the local mass (solid line) in a protoplanetary 

disk near its orbit. The migration timescale increases if the ratio 
B = A-kTioR?IM 2 < 1, which is the case for close orbits. 



a[AU] 


Figure 12. The migration timescales for a circumbinary pro¬ 
tostellar disk, based on the typical binary properties of our sim¬ 
ulations. The unperturbed Type II migration timescale, 7-niig,o 
(dashed line) shows the timescale assuming that the secondary 
mass is much smaller than the local disk mass. A longer timescale 
Tmig,i (solid line) is expected if the secondary mass is large com¬ 
pared to the disk mass (see previous Figure). 


formation rates inferred from our simulations be taken as a 
conservative estimate, with possible additional contributions 
from channels other than stellar scattering. 


4.3 X-ray output 

As discussed in HMXBs are believed to be a major source 
of X-rays in the early universe. Observations of nearby 
star-forming galaxies suggest that their X-ray luminosities 
(which are dominated by HMXBs) scale linearly with their 


star formation rate (Grimm, Gilfanov & Sunyaev|2003 

Gil- 

fanov, Grimm & Sunyaev||2004| |Persic et al.||2004| |Mirabel 

et al.||2011|). 

Mineo, Gilfanov & Sunyaev 

(2012 

1 find many 
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The first HMXBs 13 


studies support a linear proportionality between the X-ray 
luminosity of HMXBs and the star formation rate (SFR) 


(Grimm, Gilfanov & Sunyaev|2003 

Gilfanov, Grimm & Sun- 

yaev 2004 Persic et al. 2004[ and 

Mirabel et al. 2011 ). In 


the linear regime, the X-ray luminosity of the local universe 


is given by Mineo, Gilfanov & Sunyaev (20121 


keV — 3 X 10^^ X 


SFR 


Mo yr- 


(19) 


where SFR is the star formation rate. 

How the ratio of X-ray luminosity to SFR evolves with 
redshift is a key question in evaluating the properties of 
galaxies and young stellar populations, and in linking the 
earliest galaxies (most of which should be actively forming 
stars, based on their mass accretion rates) with their X-ray 
luminosities. The question is as yet unresolved by current 


observations (Dijkstra et al. 2012 Basu-Zych et al. 20131, 


and is often treated as a free parameter in studies estimating 


the X-ray production of the first galaxies (e.g. Furlanetto 
2006[ IFialkov, Barkana fc Visbal||2014| [Tanaka, O’Leary fc 


Perna||2015|)T 

In the following, based on the results of our simulations, 
we are going to quantitatively evaluate the relation between 
X-ray luminosity and SFR, and compare it with equation 
(191. We can write the X-ray luminosity as 


f'2-io keV — Tfidd X /sdd X /2-IO keV 

X face X fsur X fesc xFhmxbxSFR. (20) 


Below, we discuss each quantity in equation (20l. 


(i) F/Edd, the Eddington luminosity, which scales with 
the typical mass of the BH engine Mbh as 1.3 x 
10^®(Mbh/ M©) erg s“\ 

(ii) /Edd, the typical ratio of the total radiative power 
emitted by HMXBs (the bolometric luminosity) to Lsdd. If 
the typical luminosity of a HMXB during an active phase 
is Eddington, then /sdd is effectively the mean duty cycle. 
Other studies have adopted values ranging from 0.1 to 0.5 
([Mirabel et al.|2011 Belczynski et al.|2008 Salvaterra et al. 


2012b I; we take as our fiducial value /Edd = 0.1. 


(hi) /2-10 keV, the fraction of the bolometric luminosity 
that is emitted between 2 and 10 keV. Observational esti¬ 


mates vary between 0.1 and 0.8 (e.g. Sipior, Eracleous & 
Sigurdsson||2003 Migliari fc Fender||20(36 1. Because the BH 

masses in Pop III HMXBs are expected to be higher than in 
present-day populations, and the peak energies of accretion 
disks scale with the mass of the central engine as 
their characteristic spectra could be somewhat softer. We 
do not expect this to significantly affect our estimates here. 

(iv) tacc, the time that a massive binary spends as a 
HMXB, with the less massive star donating mass to the 
more massive BH companion. If the two stars form simulta¬ 
neously and form a compact binary before the more massive 
member dies to become a BH, then this is simply t2 — ti, the 
difference in the lifetimes of the stars. We use this (some¬ 
what arbitrary) estimate. However, because the lifetimes of 
the less massive star {t 2 > 10 Myr for stars with masses 
< 10 M©) are comparable to the expected specific star for¬ 
mation rate in galaxies at this redshift, we argue that any 
prescription that satisfies tacc < t2 is a reasonable order-of- 
magnitude estimate. 

(v) /sur, the fraction of HMXB candidates identified in 


our simulations that actually survive to become HMXBs. 
This quantity accounts for possible disruptions of binaries, 
due to (a) the merger of the stars during main sequence 
and post-main sequence evolution (Power et al. 2009| ); (b) 
the more massive star getting kicked following a supernova 


explosion (e.g. Repetto, Davies & Sigurdsson 2012 Janka 


20131; and (c) subsequent disruptions by stellar scatterings 


that were not captured by our simulations. (Our simulations 
follow the evolution of star groups until the formation of a 
stable compact binary, but not until stellar death millions of 
years later.) Theoretical estimates typically yield ~ 0.2 —0.3 
(e.g. Jeon et M]|2014 Artale, Tissera fc Pellizza|2015 1 

(vi) /esc, the fraction of X-rays that escape the galaxy. 
Unless the environment of the HMXBs are Compton-thick, 
which is unlikely for the low-mass galaxies of interest, we 
expect /esc < 1. 

(vii) Fhmxb, the number of HMXBs formed per stellar 
mass. This is the main output of our simulations. Whereas 
previous theoretical works had arrived at this value by ex¬ 
trapolating the locally observed value with an assumed red¬ 
shift evolution, or with free parameters, here we provide an 
estimate based on suites of A^-body simulations whose ini¬ 
tial conditions are motivated by cosmological simulations of 
Pop HI star formation. (Note that this quantity has units 
Mg^; we use the capital letter to distinguish it from the 
dimensionless fractions represented by /) 

Across all of our small-scale simulations—varying the 
number of stars in the group, whether groups evolved in 
isolation or through several different orientations of mu¬ 
tual collisions, and exploring two values for the ambient 
gas density separated by two orders of magnitude—we find 
Fhmxb ~ 10“^, varying by less than a factor of 4 between 
the lowest and the highest values (see Table [^. 

Finally, we can write L 2-10 keV as follows : 


F/ 2-10 keV _ r, no w -^BH /Edd /2-IO keV face 

SFR M© 0.1 0.1 Myr 

fsuT fesc Fhmxb 

^ 05 ^ OA ^ 10-3 M/i 




10^^ erg s~^ 
M© yr-i 


( 21 ) 


Based on our choices of the factors, /Edd = 0.1, /2-10 kev = 
0.1, /sur = 0.5 and /esc = 0.5, we can estimate the normal¬ 
ized X-ray luminosities per SFR, L2-iokev/SFR. We report 
this quantity for each of our models in Table|^ It varies from 
a minimum of 37 to a maximum of 450 among the studied 
scenarios. 

These Lx-to-SFR ratios are ~ 40—150 higher than what 
is observed in the local Universe. This result is qualitatively 
consistent with the findings of [Basu-Zych et al.| ( |2013| ) and 
Kaaret (20141, who find an increase in the Lx-to-SFR ra¬ 


tios toward 2^4. Our high Lx-to-SFR values stem from 
the large mass of the HMXBc primary, and the relatively 
low mass of the secondary. The former leads to a higher 
Eddington luminosity compared to typical stellar-mass BHs 
(~ 3 M©) in the local Universe, and the latter results in 
long stellar lifetimes, which in turn leads to longer face. Note 
that we used the stellar mass as a proxy for the BH mass for 
simplicity, due to the theoretical uncertainties in evaluating 
the mass loss due to winds and during the transition to a 
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BH. Any significant mass loss (e.g. simulations by Zhang, 


[Woosley fc Heg^|2008| suggest that the BH remnants of 
massive metal-free stars end up with ~ 40% of the progen¬ 
itor mass) would be directly translatable to the estimated 
X-ray luminosities reported in Table 


4.4 Implications for the thermal history of the 
IGM and the 21cm radiation 


A higher Lx-to-SFR ratio implies that IGM heating will oc¬ 
cur earlier than commonly thought. The thermal history of 
the IGM can be probed in the 21 cm line, which is observ¬ 
able in absorption (or in emission), depending on whether 
the spin temperature of the IGM is below (or above) the 
GMB temperature. 


If the IGM heats early, as suggested by our estimates 
of the X-ray emission of early galaxies, the 21 cm absorp¬ 
tion line appears earlier, and the “dip” as a function of red- 
shift caused by adiabatic cooling is not as deep and not as 
sharp as in the case of late heating as it would otherwise 
(see Figure 2 in Fialkov & Barkana 20141. Another conse¬ 
quence of early, intense heating is that the temperature of 
the IGM could become high enough, to suppress the for¬ 


mation of low-mass galaxies I Ripamonti, Mapelli & Zaroubi 


20081 and the growth of their nuclear BHs (Tanaka, Perna 
fc Haiman|2012 1. 


We can use our results for the formation rates of Pop III 
HMX Bs to estimate the fra ction of LGRBs from Pop III 
starj^ [Bromm fc Loeb (20061 quantified the GRB formation 
efficiency as 


^GRB — ^BH ^bin ^close ^beaming 5 


( 22 ) 


where ?7 bh is the number of BH-forming stars resulting 
from a given total stellar mass, r^bin is the binary frac¬ 
tion and J^ciose is the fraction of sufficiently close bina¬ 
ries to undergo RLOF. For Pop I/II stars they calculated 
?7bh — 1/(700 Mq). Combining this value with adopted 
values for the other parameters—r^bin ~ 0.5, ?7ciose ~ 0.3, 
and J?beammg (1/50) - (1/500)—yields r^GRB.Popi/ii ~ 
4.2 X (10~'^ - 10~'^) M~^ 

[Brornm fc Loeb| ( |2006[ ) noted that it was only to make 
educated guesses for the Pop III case, due to an absence of 
detailed calculations for the fraction of close binaries. Our 
work is a first attempt to fill the gap in our theoretical knowl¬ 
edge. We can write Fhmxb = »7bh Jfbin Jfdose. Adopting for 
comparison the same value of 77beaming — (1/50) — (1/500), 
we then infer ?7 grb,Pop hi — 4.8 x (10“® ~ 10“^) Mg^ for 
the interacting 5-star case, and 2.2 x (10“® ~ 10“®) Mg^ for 
the (most favorable) 10-body scenario. 

Therefore, our results suggest that LGRB rates from 
Pop III stars could be comparable to or somewhat higher 
than the rates from Pop I/II stars. 


4.5 Implications for Gamma Ray Bursts from 
Pop III stars 


As discussed in Q the fraction of HMXBs at high redshifts 
has potential implications for the expected rates of LGRBs 
from Pop HI stars. According to the collapsar model ( |Mac-| 
Fadyen fc Woosley||1999 1, for an exploding massive star to 


yield a GRB, several conditions need to be satisfied, namely: 


(i) The core of the star must collapse to a BH. This is 
realized by most Pop HI stars, given their large masses. 

(ii) The hydrogen envelope of the progenitor star must be 
stripped, so that a relativistic jet can penetrate and exit the 
remaining envelope. 

(hi) The BH should be surrounded by an accretion disk 
of high angular momentum material. This is realized if the 
core of the progenitor star has retained sufficient angular 
momentum during the evolution. 


Binary systems more easily satisfy the last two conditions 
with respect to single stars (see e.g. Cantiello et al.|[2007 |. 
In fact, for single stars to end their lives as LGRBs, they 
need to be born with large initial rotation (since they are 
less likely to be spun up by other stars), and also avoid being 
slowed down by magnetic torques (e.g. Spruit 2002 Yoon, 


Langer fc Norman|2006 Perna et al.||2014l. 

In contrast, binary stars can spin up the helium core 
of the progenitor star via tidal coupling and spin-orbit lock¬ 
ing. Further, RLOF can strip the hydrogen envelope dur¬ 
ing a common-envelope phase without reducing the rota¬ 
tion of the helium core ( [Brornm fc Loeb||2006[ ). This is es¬ 
pecially important for Pop HI stars, whose heavier hydro¬ 
gen envelopes would be more difficult to shed in isolation. 
Therefore, compact binary systems, or HMXBs, constitute 
a promising channel to produce LGRBs from Pop HI stars. 


4.6 Caveats 


Our suite of A^-body simulations, spanning diverse sets of 
initial conditions for Pop HI stars and their environments, 
point that HMXBs form in higher fractions in the earliest 
galaxies than at low redshifts, and make significant contribu¬ 
tions to the thermal history of the IGM, the 21 cm signature 
at 2 ~ 20, and to the rates of LGRBs. The fact that the for¬ 
mation rates of HMXBs varied little between the simulations 
suggest that our estimates for Fhmxb is reasonably robust. 
However, we here point out several uncertainties of our work 
that could affect our conclusions. 

One important factor is the IMF of the stars (e.g. Hi 
rano et n)||2014 |. We adopted the IMF of Stacy fc Brornm 
120131, which were based on the masses of protostars roughly 


5000 yr after the formation of the protostellar seeds. 

However, our key results are based on fragmentation of 
the protostellar clouds on smaller scales, as found in the sim¬ 
ulations by jGreif et al.| ( [^011| ). In those simulations, the most 
massive protostars had the highest accretion rates, suggest¬ 
ing that the IMF slope may be steeper than what we as¬ 
sumed. We also did not account for changes in the masses of 
stars as they evolved. These are important considerations. 


^ Note that in our study, since we are considering HMXB ’can¬ 
didates’, there is also the possibility of merger of the two stars 
during the common-envelope inspiral phase, when they both are 
stripped down to their Helium cores. This event could provide 
another avenue for the formation of GRBs ( [Fryer fc Heger|2005| l. 
Alternatively, after both stars have undergone the SN explosion, 
if the system is still bound, the compact objects of the binary, 
upon merger as a result of gravitational energy loss, would be 
likely contributors to the population of Short Gamma-Ray Bursts 
(e.g. Narayan, Paczynski fc Piran|1992^. 
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5-body (ng) 

5-body(n4) 

See. 1 

See. 2 

See. 3 

See. 

See. 4B 

Mi[ M 0 ](Tiife,i[ Myr]) 

110(2.5) 

88(2.9) 

110(2.5) 

120(2.4) 

110(2.4) 

160(2.0) 

110(2.4) 

M2 [Mq] (Tiife,2[Myr]) 

12(17) 

11(18) 

18(11) 

45(4.8) 

28(7.1) 

63(3.7) 

63 (3.7) 

tacc(= riife,2 “ Tufe, i) [Myr] 

14 

15 

8.6 

2.5 

4.7 

1.7 

1.3 

Fhmxb [ 10 “"* Mg^] 

4.6 

4.6 

15 

11 

9.0 

4.2 

8.4 

T 2 -I 0 kev/SFR [1039ergs-l M”! yr] 

220 

200 

450 

107 

160 

37 

40 

r / r local 

r>2-10keV/r^2-10keV 

75 

67 

150 

36 

52 

12 

13 

^GRB.PopIII 

’?GRB,PopI/II 

2.2 

2.2 

7.0 

5.4 

4.3 

2.0 

4.0 


Table 4. Summary of the results of X-ray luminosity and GRB efficiencies. In the table, second and third columns correspond to 
5-body calculations and the rest columns (from See. 1 to See. 4B) are the results from 10-body calculations with the number density 
of 10®cm“^.According to the ratio ^2-10 kev/^2°—10 keV’ ^ 2 °-io keV = 3 X 10^® erg s ^, it is expected that the 2-10 keV X-ray 

luminosity is ~ 10^ times larger than the X-ray luminosity of the local universe. These larger values essentially result from the large mass 
of the star 1. In the bottom row, we have listed the LGRB efficiency, derived from each set of simulations. The higher HMXB formation 
rates, lead to larger efficiencies than the estimated value for Pop I/II stars, i?GRB,Popl/ll — 4.2 X 10“® with beaming factor of 1/50 (e.g. 
[Brornm fc Loeb|2006| |. 


as the masses of the stars play a critical role in determining 
the time that a binary spends transferring mass as a HMXB. 

We made an effort to acconnt for this limitation by using 
free parameters snch as /aur, which accounts for the mass loss 
during the SN explosion. 

Another factor that conld be more carefnlly treated in 
future studies is the spin of the star. In a binary system, 
this impacts the circularization of the orbit, followed by 
low eccentricity and the synchronization of spin with or¬ 
bital phase[^ Due to tidal dissipation and circularization, 
this could additionally decrease the orbital distance of a bi¬ 
nary, boosting the formation of HMXBs and change the av¬ 
erage value of the eccentricity. 

In addition to reducing the orbital separation, the spin 
makes a difference in mass transfer rates. 

Our simulations generally produce binaries with large 
eccentricities. For highly eccentric orbits, mass transfer oc- 
pericenter | |Lajoie fc Sills|[2011| and |Sepin-| 
I, and the rates depend on whether the or¬ 
bital angnlar speed of the star is super-sychronons or sub- 
synchronous with rotational angular speed | |Davis, Siess &l\ 
Deschamps||2013 b 

We found that the viscons timescales in our eccentric 
binaries were longer than the orbital timescales in the ma¬ 
jority of cases, and used this fact to conclude that their 
duty cycle should be the same as for nearly circular binaries. 
However, we did find one exception, in which a RLOF ac¬ 
cretion event wonld be snfficiently short-lived to be episodic. 
A more careful study focused on stellar rotation effects may 
be necessary to more conclusively estimate the duty cycle of 
eccentric Pop HI HMXBs. 

Finally, our simulations are an attempt to model the 
stellar dynamics as a gravitational Wbody problem with 
perturbative forces due to a fixed, smooth gaseous back¬ 
ground. More detailed simulations that include detailed stel¬ 
lar feedback, as well as the dynamics and thermodynamics 


curs only near 


sky et al. 2010 


^ These effects had been often assumed to be due to the tidal in- 
terations with accreting gas, but recent studies suggest that the 
orbital semimajor axis and eccentricity can either increase or de¬ 
crease depending on the binary properties at pericenter | |Sepinsk^ 
Willems fc Kalogera|2007 and Sepinsky et al.|2009 \. 


of the gaseous environment, could lead to additional revela¬ 
tions about the early evolution of Pop HI star groups. 


5 SUMMARY 


In this study, we used Y-body simulations of the first stars 
to explore the formation, evolution, disruption and energy 
output of Pop HI HMXBs. The code includes gravitational 
scattering of stars, dynamical friction, and the gravitational 
potential of ambient gas. 

The initial conditions for the simulations (i.e. IMF, typ¬ 
ical star separation in the host haloes, ambient densities) are 
taken from two different sets of cosmological simulations of 
Pop III formation, namely by Stacy & Bromm (|2013| (’large- 
scale’, i.e. a few thousands of AU), and Greif et al. (20121, 
(’small-scale’, i.e. a few tens of AU). These provide two com¬ 
plementary sets in that they explore different physical scales 
for star formation (for details, see |2.4[ ). For each of the two 
scenarios, we investigated star evolution in two backgound 
gas densities, a high-density case (10®cm“^), and a lower- 
density one (10’^cm~^). 

Based on the handful of protostars per halo that are 
found in the works quoted above, we simulated systems with 
5 stars and systems with 10 stars. We found: 


(i) 5-body simulations-. If stars form in quasi-Keplerian 
disk configurations with initial separations of hundreds of 
AU, HMXBs are highly unlikely to form. In contrast, if 
stars form in compact groups separated by ~ 10 AU, as is 
expected from turbulent fragmentation, stellar scatterings 
lead to a significant HMXB formation rate. In particular, 
we found that HMXBs form at a rate of a few per 10"^ Mq 
of stars formed, independent of the ambient gas density. 

(ii) 10-body simulations: We simulated 10 stars on sep¬ 
arations of ~ 10 AU, and evolved them as isolated quasi- 
Keplerian disks, or as two colliding groups with 5 stars each. 
For the latter, we ran several different sets of collision ge¬ 
ometries. 

We found that the HMXB formation rate was a factor 
~ 1 —3 times higher than for the 5-body simulations, mainly 
due to the fact that the larger number of stars allowed for 
more hardening via stellar scattering. 


All of the small-scale simulations suggest an X-ray lumi- 
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nosity per unit star formation that is a factor ~ 10^ higher 
than what is observed in the local Universe (under the as¬ 
sumption that other variables such as the X-ray escape frac¬ 
tion from galaxies and the duty cycle of HMXBs do not differ 
significantly). These results are mostly due to the large mass 
of the most massive star of the HMXBc compared to that of 
the companion star, implying both a large face as well as a 
higher luminosity of the remnant BH. The fact that we found 
little variation in this quantity across all of our simulations 
suggests that this is a robust estimate. Additional factors, 
such as in-disk migration of nascent stars, could further in¬ 
crease the HMXB formation efficiency. 

A direct consequence is that X-rays can heat the IGM 
rapidly at Cosmic Dawn. Signals of early heating can be 
probed via the 21 cm line radiation: the absorption line sig¬ 
nal is expected to show a broader, shallower minimum due 
to the shorter gas cooling time, while the emission line would 
be observed earlier because of the higher gas temperature, 
at earlier times. Several studies have modeled the 21 cm sig¬ 
nature of the first HMXBs, but relied on assumptions for 
their Lx/SFR relation relative to the empirical value found 
at lower redshifts. Our work provides a theoretically driven 
estimate for this quantity. 

In addition to the implications for the thermal history 
of the IGM, these high formation rates of HMXBs per stellar 
mass imply a higher GRB formation efficiency from Pop HI 
stars in binaries. This predictions can be tested with a long 
baseline of observational data from Swift. 
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